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Improved  Fortran  Program  for  Single  Particle  Energy  Levels  and 
Wave  Functions  in  Nuclear  Structure  Calculations 

Randall  S.  Caswell 


An  improved  program  has  been  developed  for 
numerical  calculations  of  single  particle  energy 
levels  and  wave  functions  for  Woods- Saxon  or  other 
real  potential  wells.   The  program  is  designed  to 
be  used  as  a  subroutine  of  a  larger  program  for 
nuclear  structure  calculation  but  may  be  used 
separately  if  desired.   Improvements  over  the 
previous  program  include:   (1)  for  bound  states 
the  wave  function  is  calculated  by  integrating 
outward  from  zero  radius  and  inward  from  a  maximum 
radius  beyond  the  nuclear  potential  well  with  match- 
ing at  an  intermediate  radius;  (2)  a  first  order 
correction  has  been  made  to  the  starting  conditions 
at  small  radius  for  the  integration  of  the  wave  func- 
tion; (3)  a  large  step  size  is  used  until  the  calcula- 
tion has  nearly  converged  on  the  eigenvalue;  then  a 
small  step  size  is  used  to  provide  maximum  accuracy; 
(4)  the  correction  of  a  small  error  in  the  Kutta- 
Runge  integration  procedure  has  been  made.   As  a  re- 
sult of  the  changes  the  program  is  approximately  five 
times  faster  than  the  previous  version.   Checks  with 
harmonic  oscillator  potentials  give  an  accuracy  of 
.01%  in  the  energy  level  values.   Listings  of  Fortran  II 
and  Fortran  IV  versions  are  given. 

KeyWords:  Eigenvalues,  eigenfunctions ,  energy  levels  , 
wave  functions,  single  particle  shell  model. 

1.   Introduction 
Nuclear  shell  model  calculations  usually  have  employed  harmonic 
oscillator  wave  functions  for  the  evaluation  of  radial  integrals  for 
both  bound  and  continuum  states.   Deeply  bound  states  can  be  well 
represented  by  harmonic  oscillator  wave  functions.   For  continuum 
wave  functions  it  appears  desirable  to  use  more  realistic  wave  func- 
tions, for  example,  those  of  a  real  Woods- Saxon  well.   A  previous 
program  (Caswell,  1962)  was  developed  to  provide  these  eigenvalues 
and  wave  functions.   The  older  program,  however,  used  a  method  of 
matching . inside  and  outside  wave  functions  at  successively  larger  radii 
which  is  inherently  slow.   In  some  large  nuclear  structure  programs 


the  calculation  of  the  single  particle  wave  functions  occupies  an 
undesirably  large  fraction  of  the  computer  time.   The  program  described 
here  speeds  up  the  calculation  by  a  factor  of  roughly  5»  provides 
increased  accuracy,  and  gives  normalized  wave  functions  at  up  to  201 
radial  points. 


2.  Description  of  the  Calculation 

The  chief  calculational  steps  are:   (l)  calculation  of  the  Woods- 
Saxon  real  potential  well  including  a  spin-orbit  term;  (2)  numerical 
integration  of  the  wave  function,  u(r)  =  r|(r),  from  zero  to  a  match- 
ing radius,  RMATCH  and,  for  bound  states,  from  a  maximum  radius  beyond 
the  nuclear  potential,  RMAX,  backwards  in  radius  to  the  same  matching 
radius  --  for  scattering  states  the  maximum  radius  is  used  as  the 
matching  radius;  and  (3)  calculation  of  the  logarithmic  derivatives, 
f  (inside)  and  f  (outside).  For  scattering  states  f  (outside)  is 
determined  from  the  criterion  of  a  90°  phase  shift.   It  is  also  possible 
to  read  in  from  the  input  data  values  for  f  .(outside)  which  are  taken 

JC 

as  the  value  at  the  maximum  radius,  RMAX„  The  approximate  location  of 
the  eigenvalues  may  be  determined  by  a  survey  versus  energy  of  f  (in- 

aj 

side),  f  (outside) ,  and  the  difference  between  them,  called  "signature 
of  the  well".  An  automatic  search  procedure  finds  the  energy  for  which 
f  (inside)  =  f  (outside),  and  the  wave  function  vs.  radius  is  calculated 

Xj  AJ 

and  stored  in  memory. 

The  potential  is  calculated  by  subroutine  VR3  as  a  function  of 
radius  at  up  to  399  points,  which  is  the  number  required  for  the  numer- 
ical integration  routine  for  up  to  201  radial  points  for  the  wave 
functions.   Spacing  between  points  is  RMAX/NINTVL  where  RMAX  is  the 
maximum  radius  and  NINTVL  is  the  number  of  intervals  (one  less  than 
the  number  of  points).   The  potential  used  is: 

V(r)  =  Vcp(r)  -  „Vc  (■£)*    <*3)i  g  (1) 

where  p(r)  =  - —  ,  the  "Saxon"  potential;  a  is  the  strength 

1  +  expf 


of  the  spin-orbit  interaction  compared  to  the  "Thomas"  interaction 

for  a  nucleon;  M  is  the  nucleon  mass;  V  is  the  central  real  potential; 

ft/Mc  is  the  Compton  wave  length  for  a  nucleon;  Z   is  the  orbital  angular 

— > 

momentum  of  the  neutron;  a   is  the  Pauli  spin  operator  of  the  incident 

neutron;  R  is  the  nuclear  radius;  and  a  is  the  "diffuseness"  para- 
o  r 

meter  of  the  potential „   If  a  negative  value  is  given  for  a>  it  is 
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taken  as  V  ,  the  depth  of  the  spin-orbit  well  in  MeV  where  a,V  (ttt" 

SO        /  a  \2  C\cdrlC 

is  replaced  by  V  ( )  ,  m  being  the  mass  of  the  jt-meson.   If 

r       J      soVm  c/    at 
it 

j  =  Z  +  \t    then  a»i  =  I,   and  the  values  of  the  variable  VPLUS  (nuclear 
potential)  for  different  radii  are  calculated.   If  j  =  JL  -  \t    then 
a*  Z  =   -(j6+1),  and  potential  is  calculated  for  the  spin-orbit  anti- 
parallel  case.  Other  potential  shapes  can  be  used  by  changes  in  this 
subroutine.  For  example,  the  program  has  been  tested  with  a  harmonic 
oscillator  potential. 

The  calculation  of  f  (inside)  and  f  (outside),  and  the  search 
procedure  for  matching  the  f  values  is  controlled  by  EIGENV,  the  main 
subroutine  of  the  program.  First,  the  radial  part  of  the  non-relativ- 
istic  Schrb'dinger  wave  equation  is  integrated  out  in  radius  for  the 
appropriate  Z: 


d2u.(r)   2m  r    Z(Z+l)ft 
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-1—  +- 

dr3     tf 


E V(r)lu  (r)  =  0.  (2) 

2mr3         J  l 


The  values  of  the  function  and  the  first  derivative  at  the  first 

radial  point  (radius  =  H  )  are  determined  from  the  relation  u  (r)  = 

JL+l 
C  r    as  r  ->  0,  where  C  is  a  constant.  We  make  a  first  order  correc- 

tion  to  the  starting  value  of  the  wave  function,  u,  and  the  deriva- 

du     ,   -. .       , .  -, 
tive  j—  at  the  first  radial  point  a; 

rewritten  (omitting  subscript  i,'s): 


tive  j—  at  the  first  radial  point  as  follows:  Equation  (2)  may  be 


£u  =  Ki+lju  m  2h  (e.v)u>  (3) 

dr3  r3  fi2 


As  r  -»  0,  the  last  term  on  the  right  is  negligible  compared  to  the 
first  and  our  solution  is  just  u  =  r  as  indicated  above,  taking 
C  equal  to  one.  Taking  into  account  the  second  term  and  using  our 
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approximate  solution  for  u,  we  may  find  the  derivative  at  the  first 
point  r  (=  HI)  by  integration  of  (3). 

du   , „  -v  I       2m  (E-V)r  m 

d7=  u+1)r  -W   (ilL)        •  w 

Integrating  again, 


_  i+1  .  2m   (E-V)     1+25  f5x 

u_r    W  (l+2)U+3)  r   *  t5; 


Equations  (1^)  and  (5)  are  used  as  the  initial  values  of  wave  function 
and  derivative  to  be  used  by  the  wave  function  integration  subroutine, 
V  is  taken  as  the  value  of  the  potential  at  the  point  r  =  H  . 

A  similar  problem  arises  in  the  backward  integration  from  a 
radius  EMAX  for  bound  states.  At  this  distance  the  centrifugal  poten- 
tial term  is  small  and  the  differential  equation  may  be  written: 

0  =  f1[V(r)-E]u  =  72u  (6) 

where  the  expression  in  the  square  bracket  is  positive  for  a  bound 
state.   If  we  consider  V(r)  as  constant,  then  the  descending  ex- 
ponential solution  to  this  equation  is  u  =  exp(-7r),  du/dr  = 
-y   exp(-7r),  and  (du/dr) /r  =  -7.  Since  normalization  is  initially 
arbitrary,  we  set  a  =  1  and  (du/dr)  =  -7  at  RMAX  for  starting  the 
backward  integration. 

The  second  order  Runge-Kutta  method  is  used  for  the  numerical 
integration  (see  Zurmuhl,  I96I) .  An  error  in  the  previous  code 
(Caswell,  I962)  which  led  to  errors  of  about  1  percent  in  the  eigen- 
values and  eigenfunctions  by  causing  a  slightly  wrong  potential  to  be 


ased  at  certain  steps  has  been  corrected.   The  integration  is  con- 
trolled by  subroutine  INTEG3  and  carried  out  by  subroutines  KRR1 
and  FR. 

Modes  of  the  calculation.   The  calculation  may  be  run  in  three  al- 
ternative modes.   When  MODE  =  1,  if  the  energy  is  negative  (bound  state), 
the  wave  function  is  integrated  outward  from  zero  to  RMATCH  and  inward 
from  RMAX  to  RMATCH,  using  the  starting  conditions  given  above.   For  con- 
vergence on  an  eigenvalue  it  is  required  that  the  "inside"  and  "outside" 
f  's  match  at  RMATCH.   If  the  energy  is  positive  (continuum  state),  the 
criterion  of  a  90  phase  shift  is  applied  to  the  wave  function  at  the 
radius  RMAX  where  the  inside  and  outside  wave  functions  must  join  smooth- 
ly (as  indicated  by  having  the  same  value  of  f.)«   For  the  90  phase 

Xj 

shift  continuum  case,  the  "outside"  or  asymptotic  solutions  are  spheri- 
cal Neumann  functions,  G.  =  -xn  (x) ,  since  Coulomb  effects  are  not 

Xj  Xj 

considered.   The  argument  x  =  kr.   The  evaluation  of  the  logarithmic 
derivatives,  f  ,  is  given  in  Appendix  1. 

Xj 

When  MODE  =  2,  values  of  f  (outside)  at  a  radius  beyond  the  nuclear 

Xj 

I  potential  radius,  RMAX,  may  be  given  the  EIGENV  program  by  the  calling 
I  program.  The  variable  FLOUT  is  used  to  read  in  the  value  of  f  (out- 

Xj 

side).   In  the  case  of  a  bound  state,  after  integration  inward  from 
RMAX  starting  at  the  value  FLOUT,  the  variable  FLOUT  is  actually 
i  matched  to  the  f  (inside). 

When  MODE  =  3,  the  program  just  calculates  and  normalizes  the  wave 
function  at  the  given  energy  EN.   No  search  is  carried  out.   EN  must 
have  been  determined  on  a  previous  run.   This  mode  is  provided  to  mini- 
mize running  time  by  avoiding  unnecessary  searches. 

3.   Search  Procedure 
To  begin  the  search  procedure,  one  needs  a  trial  energy.   In 
principle  this  may  be  selected  by  guess.   The  guess  will  usually  con- 
verge to  some  eigenvalue  if  enough  search  steps  are  allowed.   A  more 
systematic  way  to  find  approximate  energies  for  starting  the  automatic 
search  procedure  is  to  use  the  method  called  the  "signature  of  the  well." 
This  is  carried  out  by  setting  the  variable  NOMAX  =  1,  and  making  a  run 
with  a  large  number  of  trial  energies,  for  example,  at  1  MeV  steps 

5 


throughout  the  energy  range  of  interest.   When  NOMAX  =  1,  the  radius  of 
the  well  at  the  energy,  EN,  is  used  as  the  matching  radius  for  bound 
states,  and  RMAX  is  used  for  scattering  states.   Values  are  calculated 
for  f  (inside),  f  (outside),  and  for  the  difference  between  the  two 

Xj  Xj 

(see,  for  example,  Figure  1).   The  integration  in  radius  is  carried  out 
only  once  for  each  value  of  EN.   The  energy  values  at  which  the  dif- 
ference, DIF,  equals  zero  are  close  to  the  energy  of  the  final  converged 
eigenvalue  which  is  to  be  found,  and  they  serve  as  suitable  starting 
values  for  the  automatic  search  procedure. 

Automatic  search  procedure.   First  an  integration  of  the  wave  func- 
tion is  carried  out  at  the  trial  eigenvalue  energy  (the  value  of  EN 
when  EIGENV  is  called).   A  second  integration  is  carried  out  at  a  nearby 
energy  (one  percent  difference  in  energy).   Using  the  differences  (be- 
tween the  inside  and  outside  f  * s  found  at  the  matching  radius)  found 

Xj 

in  these  two  trial  calculations,  a  linear  interpolation  or  extrapolation 
is  made  to  an  energy  for  which  the  difference  is  predicted  to  be  zero 
(corresponding  to  smooth  joining  of  the  inside  and  outside  wave  func- 
tions).  Integration  in  radius  is  carried  out  at  this  predicted  eigen- 
value energy,  and  the  difference  is  again  found.   We  now  know  three 
energies  and  the  three  corresponding  differences.   Using  this  informa- 
tion, the  program  predicts  by  quadratic  interpolation  a  fourth  energy 
at  which  the  difference  should  be  closer  to  zero.   The  calculation 
continues,  energy  prediction  being  made  by  quadratic  interpolation  on 
the  last  three  values,  until  either  the  difference  between  f.'s  is 

Xj 

closer  to  zero  than  the  difference  specified  by  the  value  of  EPS  given 
to  the  subroutine,  or  until  the  maximum  number  of  tries,  NOMAX,  has 
been  used  up.   The  wave  functions  during  the  convergence  procedure  are 
shown  in  Figure  2  for  a  typical  case.   If  the  maximum  number  of  tries 
has  been  used  up  without  converging  on  an  eigenvalue,  the  value  of  EN 
and  values  of  the  wave  function  at  radial  points,  WAVEFN(I),  are  set 
equal  to  zero  to  insure  that  the  calling  program  does  not  use  incorrect 
information.   When  the  converged  eigenvalue  has  been  found,  the  part  of 
the  wave  function  between  RMATCH  and  RMAX  is  renormalized  so  that  it  has 
the  proper  absolute  magnitude  to  join  smoothly  to  the  inner  wave 
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function  at  RMATCH.   A  radial  integration  of  the  Ju3dr  is  then  carried 
out  by  subroutine  SIMPSN,  and  the  complete  wave  function  is  normalized 


so  that 


subroutine  ill 
J  u2dr  =  1. 


o 
To  save  computer  time  an  integration  step  size  five  times  the  final 

step  size  is  used  for  the  initial  steps  of  the  search  (when  DIFX5.1). 

Integration  step  size  is  controlled  by  the  variable  KSTEP  which  is  set 

jequal  to  5  in  the  early,  rough  part  of  the  search  and  set  equal  to  1  for 

the  final  determination  of  the  eigenvalue  and  wave  function. 

4.   Sample  Results 

Sample  results  for  eigenvalues  and  wave  functions  are  shown  in 
Figure  3.   Eigenvalues  are  shown  for  the  lsi,  2si  ,  3si,  lp^»  lPi»  Ids 
state.   Wave  functions  are  shown  for  the  15i,  25i,  and  Ids.  states. 
5.   Accuracy  and  Checks  of  the  Code 

To  test  the  accuracy  of  the  program,  subroutine  VR3  was  rewritten 
to  use  an  harmonic  oscillator  potential: 


V(r)  =  -100  I  1  -  (r/5)3 


MeV 


where  r  is  in  fermis.   For  this  particular  potential  hco  =  18.810  MeV 
Ifor  AMASS  =  15.0.  Based  on  this  value  we  may  compare  eigenvalues  given 
by  the  code  for  various  I   values  with  the  usual  relation  E  =  (n+j-Ohoo, 
all  energies  being  in  MeV. 

Harmonic  Oscillator 

Levels s  states        p  states        d  states 

-15.355  -15.352(2p) 

-34.165  -34.162(2s)  -34.l63(ld) 

-52.975  -52.974(lp) 

-71.785  -71.784(ls) 

It  can  be  seen  that  the  code  predicts  the  level  separations  to  an 
accuracy  of  about  0.01%  which  is  about  the  accuracy  of  various  constants 
in  the  code. 
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Appendix  1.      Evaluation  of  f    (outside)    for  Scattering  States 

For  continuum  states    (EX))  we  may  write  for  the   logarithmic   deriva- 


tive: 


(du  /dr)              [~(d/dr)(rn.)lR 
f^(outside)   =  R £ ^  =  R  ±= Ui        ,  (!) 

(uA  [rnA 


where  n,  is  the  spherical  Neumann  function,  which  corresponds  to  a  90° 

phase  shift,  as  defined  in  Schiff  (1955). 

But 


—  (xn^)  =  —  (x*+1x~^(x))  =  (jH-l)[nA(x)]  -  **mW.      (2) 


/2m  xi 
If  X  =  kR  =  (— ^  E)  R,  we  have  for  f , (outside) 

V  J  l 


Xn4^(X)  ,  6/  . 

(3) 


f,  =  jH-1 ^± (=  x-^-1, 

I  Vx>  G/ 

where  the  prime  refers  to  differentiation  with  respect  to  X.   The 
quantities  F  and  G  are  the  "regular"  and  "irregular"  wave  functions 
which  for  neutrons  (no  Coulomb  interaction)  are  F  =  X  j.(X)  and 

Gi =  -x  vx)- 

For  values  of  X>1,  "forward"  recursion  formulae  are  used  (see, 
for  example,  Abramowitz  and  Stegun,  1964): 

cosX 

nQ(X)  = , 

X 
cosX   sinX 

n,(X)  = 5 ,  and 

XX 

n  (X) 
nx+1(X)  =  (2X+1)-—  ■  nM(X). 


For  I  <■   9  the  values  accurately  reproduce  tabular  values  from  the  same 
reference,  there  appearing  to  be  no  need  for  use  of  "backward"  recursion 
formulas  (from  large  to  small  I). 

For  X  <  1  a  series  expansion  is  used  (see  Abramowitz  and  Stegun, 
1964) : 


1-3' 5'...  (2j&-1)  ,     h?  (|x2)a 

n/X)  =        "1+1         \}~  11(1-21)    +  2l(l-2JL)(3-2Z)   "  "  J  * 
X 
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Appendix  2.   Instructions  for  use  of  EIGENV  subroutine  and  meaning  of 

input  variables. 

To  call  the  EIGENV  subroutine  values  must  be  assigned  to  the  follow- 
ing list  of  variables  (with  meanings  given)  and  then  the  statement  CALL 
EIGENV  calls  the  subroutine.   The  subroutine  returns  as  the  answer  a 
value  of  EN  which  is  the  energy  eigenvalue,  and  of  a  normalized  WAVEFN(I) 
at  up  to  201  equally  spaced  points.   The  wave  function  is  rf. 

TYPICAL  VALUE   MEANING 


3.15 


radius  of  the  Saxon  potential  in  fermis, 
typically  1.25  A1/3. 


.65 


diffuseness  parameter  of  the  Saxon  poten- 
tial in  fermis. 


50.0 


central  real  potential  in  MeV  (is  negative 
for  attractive  potentials). 


20.0 


magnitude  of  the  spin-orbit  potential  ex- 
pressed in  number  of  times  larger  than  for 
the  Thomas  term  of  a  nucleon.   If  a  negative 
number  is  given,  it  will  be  taken  as  V   the 
depth  of  the  spin-orbit  terra  in  MeV. 


so 


15.0 


atomic  mass  in  a.m.u.  of  the  nucleus  exclud- 
ing the  single  neutron  whose  states  are 
being  considered. 


8.0 


maximum  radius  to  which  wave  function  will 
be  calculated  in  fermis  (typically  about 
5  fra  greater  than  the  average  radius  of  the 
well). 


,0001 


the  difference  between  the  value  of  f .  (in- 
side) and  f  .(outside)  must  be  less  than  EPS 
for  the  eigenvalue  to  be  considered  deter- 
mined. 


normal  search  for  bound  or  scattering  states, 
If  a  bound  state,  f  .(outside)  is  determined 
by  integration  inward  from  RMAX  to  the 
matching  radius  RMATCH  (which  is  automati- 
cally determined).   For  a  scattering  state 
the  matching  is  at  RMAX  and  f  .(outside)  is 
determined  from  the  criterion  for  a  90° 
phase  shift. 


f  .(outside)  is  read  in  from  the  input  data. 
It  is  considered  as  the  value  of  f  .(outside) 
at  radius  RMAX. 
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VARIABLE 


MODE 


NOMAX 


TYPICAL  VALUE   MEANING 
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just  calculate  the  wave  function  for  the 
given  energy,  EN.   Do  not  search  for  a 
better  eigenvalue. 

maximum  number  of  automatic  search  steps 
seeking  the  eigenvalue.   If  N0MAX=1,  one 
point  is  determined  for  the  "signature  of 
the  well." 


■ 


NINTVL 


100 


JDBLD 


JPRINT 


EN 


1 
2 

11.0 


number  of  intervals  in  radius  over  which  the 
wave  function  is  to  be  determined.   In  the 
example  101  radial  points  will  be  calcula- 
ted, equally  spaced  between  the  first  point 
at  R=0  and  the  last  point  at  R=RMAX. 
NINTVL  must  be  a  multiple  of  10. 

orbital  angular  momentum  of  the  single 
particle  state  sought. 

twice  the  j -value  of  the  state  sought.   In 
the  example  j=-f  so  "j-doubled"  is  3. 

no  printing  of  intermediate  results. 

print  intermediate  results  for  code  check- 
ing. 

value  given  to  subroutine  is  the  initial 
energy.   Value  returned  by  subroutine  is 
the  eigenvalue  (if  successful). 
Value  is  in  MeV. 


FLOUT 


4. 


read-in  value  of  f  .(outside)  for  the  radius 

RMAX. 
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Appendix  3.   Fortran  II  listing. 
The  subroutine  and  their  functions  are: 

MAINEV  A  main  program  for  testing  the  EIGENV 

program  and  a  sample  for  how  to  call  the 
program. 

EIGENV  Controlling  subroutine  which  performs  match- 

ing and  eigenvalue  prediction. 

VR3  Calculates  V(r). 

INTEG3  Controls  integration  in  radius. 

KRR1  Kutta-Runge  integration  subroutine. 

FR  Calculates  f(r)  for  differential  equation 

d3u/dra  =  f(r). 

SIMPSON  Simpson  rule  integration  subroutine  for 

normalization  of  the  wave  function. 

The  DIMENSION  and  COMMON  statements  appearing  in  the  subroutines 
must  also  appear  in  the  calling  program. 
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*      LABEL 

CMAINEV 

»      LIST 

C      MAIN  PROGRAM  TO  TEST  EIGENV  SINGLE  PARTICLE  EIGENVALUE  AND 

C      WAVE  FUNCTION  PROGRAM 

DIMENSION  RPRINT(5) , WAVE (5) ,ENRGY( 10O),FL(lO0) 
DIMENSION  WAVEFN(201),VPLUS(3  99),XINTGD(201) 

COMMON  WAVEFN, RADIUS, A, VC, ALPHA, AMASS, RMAX, EPS, MODE, NOMA X,NI NT VL  , 
1L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX, 
2F,RMATCH,I , XPLUS , XPLUSP ,FRR, INOUT, KSTEP , REFL , XINTGD , YINTG,XI NT  1 ,H, 
3XINT2,X 
5  WRITE  OUTPUT  TAPE  6,7 
7  FORMAT(lHl) 

READ  INPUT  TAPE  5,10 
10  F0RMATI72H 

1  ) 

WRITE  OUTPUT  TAPE  6,10 
WRITE  OUTPUT  TAPE  6,4501 
4501  F0RMAT(70H     RADIUS  A         VC   ALPHA(VSO)    AMASS 

1RMAX        EPS  ) 
READ  INPUT  TAPE  5 , 1 , RAD IUS, A, VC , ALPHA  , AM ASS , RMAX, EPS 
WRITE  OUTPUT  TAPE  6 , 1 ,RAD IUS, A, VC, ALPHA , AMASS, RMAX, EPS 
1  F0RMAT(7F10.5) 

WRITE  OUTPUT  TAPE  6,4503 

4503  F0RMAT(70H   ENERGIES       MODE      NOMAX     NINTVL  L      J 
1DBLD     JPRINT  ) 

READ  INPUT  TAPE  5 ,20 , KZ , MODE, NOMAX, N INT VL ,L , JDBLD , JPRINT 
WRITEOUTPUTTAPE  6 ,2  0, KZ , MODE, NOMAX, M NT VL , L , JDBLD, JPRINT 
20  F0RMAT(7I10) 

NP0INT=NINTVL+1 

GO  TO  (4081, 4082, 4081), MODE 

4081  WRITE  OUTPUT  TAPE  6,4504 

4504  F0RMATU5H  TRIAL  ENERGIES  ) 

READ  INPUT  TAPE  5 , 1 , ( ENRGY ( KY ) , KY=1 ,KZ ) 
WRITEOUTPUT  TAPE6 , 1 ,( ENRGY < KY ), KY=1 ,KZ) 
GO  TO  4085 

4082  WRITE  OUTPUT  TAPE  6,4083 

4083  F0RMAT(25H   ENERGY        FLOUT, RMAX  ) 

READ  INPUT  TAPE  5,    1 ,< ENRGY <KY) , FL ( KY ) ,KY=1 ,KZ ) 

WRITE  OUTPUT  TAPE  6,    1, ( ENRGY ( KY ) ,FL ( KY ), KY=1 ,KZ ) 
4085  DO  100  KY=1,KZ 

EN=ENRGY(KY) 

FLOUT=FL(KY) 

WRITE  OUTPUT  TAPE  6,6180 
6180  FORMAT(  68H  ENERGY     DIFFERENCE      FL  INSIDE       FL  OUT 

1SIDE    RMATCH  ) 

CALL  EIGENV 

IF(ABSF(EN)-1.E-10)  100,4088,4088 

4088  WRITE  OUTPUTTAPE  6, 4089, EN 

4089  FORMAT  (34H   WAVE  FUNCTION  VS  RADIUS,  ENERGY=, F10. 5 , 4H  MEV) 
R=-HI 
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INMAX=5 

DO  200  N=1,NP0INT,5 

DO  150  IN=1,INMAX 

NN=N+IN-1 

WAVE(IN)=WAVEFN(NN) 

R=R+HI 

RPRINT{  IN)  =R 

IF(NN-NPOINT) 150,150,155 

150 

CONTINUE 
GO  TO  159 

155 

INMAX=IN-1 

159 

WRITE  OUTPUT  TAPE  6 , 160, ( RPR  I  NT ( IN ) , WAVE ( IN ) , 

IN  =  1 

., INMAX) 

160 

F0RMAT(5(F9.3,E15.6)) 

200 

CONTINUE 

GO  TO  (100,170),JPRINT 

170 

WRITE  OUTPUT  TAPE  6,180,YINTG 

180 

F0RMAT(7H  YINTG= , F10.5 ) 

100 

CONTINUE 

READ  INPUT  TAPE  5,20,NDUMMY 

IF  (NDUMMY) 1000, 1000,5 

1000 

CALL  EXIT 

END 

• 

LABEL 

CEIGENV 

*      LIST 

SUBROUTINE  EIGENV 
DIMENSION  Y(ll) 

DIMENSION  WAVEFN(201),VPLUS(399),XINTGD(201) 

COMMON  W AVEFN, RADIUS, A, VC, ALP HA, AMASS, RMAX, EPS, MODE, NOMAX,NI  NTVL  , 
1L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX, 
2F,RMATCH,I,XPLUS,XPLUSP,FRR,IN0UT,KSTEP,REFL,XINTGD,YINTG,XINT1,H, 
3XINT2,X 
K=L  +  1 

HI=RMAX/FLOATF(NINTVL) 
NP0INT=NINTVL+1 
IF(ALPHA)  1,2,2 

1  ALPHA=ALPHA*181.1/VC 

2  CALL  VR3 
KSTEP=1 
H=HI 
GO  TO  <15, 15, 6130), MODE 

15  IF(EN)  180,180,170 
170  RMATCH=RMAX 

GO  TO  190 
180  RMATCH=(5.*FL0ATF(NINTVL/15+1))*HI 
190  MPRINT=2 

IF(NOMAX-l)  25,30,25 
25  KSTEP=5 

H=5«*HI 
30  CALL  INTEG3 

N0  =  1 
203  ENERGY=EN 

IF<EN)  213,213,220 
213  FLOMT=R«XPLUSP/XPLUS 

DIF=REFL-FLOMT 

GO  TO  6001 
220  PM0M=SQRTF(.04826*F*EN) 

DIF=REFL-FL0UT 

FL0MT=FL0UT 

GO  TO  (2201, 6001, 6130), MODE 
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2201  XX=R*PMQM 
KPLUS1=K+1 
IF(XX-1.)  300,300,400 

300  NPASS=1 

301  SUM=1. 
TERM=1. 

DO  310  N=l,20 

TERM=-TERM»0.5*XX**2/(FL0ATFtN)»FL0ATF(2*N-l-2*L) ) 

SUM=SUM+TERM 

IF(ABSF(TERM)-1.E-8)320,310,310 
310  CONTINUE 
320  IF(K-l)  325,335,325 
325  00  330  IK=2,K 

AK=IK 
330  SUM=SUM*(2.*AK-3„) 
335  GO  Tl)  (340,350),NPASS 
340  YL=-SUM/(XX**K) 

K=K  +  1 

L  =  L  +  1 

NPASS=2 

GO  TO  301 
350  YLPLUS=-SUM/(XX**K) 

K=K-1 

L  =  L-1 

GO  TO  500 
400  Y( 1)=-C0SF(XX)/XX 

Y(2)=-C0SF(XX)/XX**2-SINF(XX)/XX 

DO  410  IK=3,KPLUS1 

AK  =  IK 
410  Y(IK)=(2.*AK-3.)*Y( IK-1 ) /XX-Y ( IK-2 ) 

YL=Y(K) 

YLPLUS=Y(KPLUS1) 
500  FL0MT=FL0ATF(L)+1.-XX*YLPLUS/YL 

DIF=REFL-FLOMT 
6001  IF(NO-NOMAX)  6004,6100,6100 
6C04  IF(N0-2)  6005,6080,6090 
6005  EN=0o99*EN 

GO  TO  6100 
6080  EN     =(ENERG1*DIF-ENERGY*DIFM1)/(DIF-DIFM1 J 

GO  TO  6100 
6090  EN     =(((ENERG2*DIFM1-ENERG1«DIFM2)»DIF/(DIFM1-DIFM2))-((ENERG1 
1*DIF-ENERGY*DIFM1)*DIFM2/(DIF-DIFMD)  )  /  (  DIF-DIFM2  ) 

6100  IF(ABSF((EN    / ENERGY )- 1. )-. 2 )  6110,6110,6101 

6101  EN= ENERGY+ ABSF ( ENERGY )*S I GNF( 0,1, (EN- ENERGY) ) 
6110  WRITE  OUTPUT  TAPE  6 , 6160, ENERGY ,DIF ,REFL , FLOMT ,R 
6160  F0RMAT(F15.5,3E15.7,F8.3) 

IF  (NO-NOMAX)  6114,6140,6140 

6114  IF(ABSF(DIF    )-EPS  )  6115,6115,6120 

6115  IF(KSTEP-l)  105,6116,105 

6116  EN=ENERGY 
GO  TO  6130 

6120  IF(KSTEP-l)  100,110,100 
100  IF(ABSF(DIF    )-.l)  105,110,110 
105  KSTEP=1 

H=HI 
110  CALL  INTEG3 

DIFM2=DIFM1 

DIFM1=DIF 

ENERG2=ENERG1 

ENERG1=ENERGY 
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N0=N0+1 

GO  TO  203 

6130 

MPRINT=1 

CALL  INTEG3 

GO  TO  6150 

6140 

EN  =  0. 

DO  6145  N*1,NP0INT 

6145 

WAVEFN(N)=0. 

6150 

CONTINUE 

RETURN 

END 

• 

LABEL 

CVR 

* 

LIST 

SUBROUTINE  VR3 

DIMENSION  WAVEFN(201),VPLUS(399),XINTGD(201) 

COMMON  WAV EFN, RADIUS , A, VCr ALPHA, AMASS, RMAX, EPS, MODE, NOMAX ,N INTVL , 
1L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX, 
2F,RMATCH,I,XPLUS,XPLUSP,FRR,IN0UT,KSTEP,REFL,XINTGD,YINTG,XINT1,H, 
3XINT2,X 

R  =  HI 

H=HI/2. 

GO  TO  (2090, 2000), JPRINT 
2000  WRITE  OUTPUT  TAPE  6,2091 
2091  F0RMAT(25H  V  R  ) 

2090  C=  ,0110270/A 

NMAX=2*NINTVL-1 

DO  2210  I=1,NMAX 

REXP=EXPF( (R-RADIUS)/A) 

VCC=VC/(1.0+REXP) 

VS=ALPHA*VC*C*REXP/( ( ( 1 .O+REXP ) **2 ) *R ) 

AK=K 

IF(JD8LD-2*L)  22,22,21 
21  VPLUSl I )=VCC+VS*(AK-1.Q) 

GO  TO  25 

122  VPLUS  (  I  )=VCC-VS*AK 
25  GO  TO  (40, 28), JPRINT 
28  WRITE  OUTPUT  TAPE  6 , 30, VPLUSt I ) , R 
30  FORMAT(  E15.7*F7.3) 
40  R=R+H 
2210  CONTINUE 
RETURN 
END 
'•      LABEL 
CINTEG 
'•      LIST 

SUBROUTINE  INTEG3 

DIMENSION  WAVEFN(201),VPLUS(399),XINTGD(201) 

COMMON  WAVEFN^RADIUS, A, VC^ALPHA, AMASS, RMAX, EPS, MODE, NOMAX, NINTVL, 
1L, ELL, K,JDBLD,JPR INT, MPRINT.NPO INT, EN, FLOUT, H I, R,RR, VPLUS, V,NMAX, 
2F,RMATCH,I,XPLUS,XPLUSP,FRR,IN0UT»KSTEP,REFL,XINTGD,YINTG,XINT1,H, 

3XINT2,X 
GO  TO  (1,15) ,MPRINT 
1  GO  TO  (5,5,15),M0DE 
5  IF(EN)  5215,4221,4221 
15  WAVEFN(1)=0. 
AK=K 

F= AMASS/ (AMASS+1. 00898) 
R=H 
RR  =  R 
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V=VPLUS(2*KSTEP-1) 

XPLUS=R**K 

X=XPLUS 

ELL=0. 

CALL  FR 

ELL=K-1 

XPLUS=XPLUS+FRR*R**2/( (  ELL+2. )*(ELL+3. 

)  ) 

WAVEFN(2)=XPLUS 

XPLUSP    =AK*(R  **(K-1))  +FRR*R/ ( ELL+2 

.  ) 

GO  TO  (4090, 4090, 4080)»M0DE 

4080 

N0MAX=2 

4090 

IN0UT=1 

DC  4210  I=KSTEP,1000,KSTEP 

1  =  1 

CALL  KRR1 

R=R+H 

4127 

IF(NGMAX-l)  4130,4128,4130 

4128 

RMATCH=RMAX 

IF  IR-0.2*RADIUS )  4130,4129,4129 

4129 

IF(EN-V)      10,10,4132 

10 

RMATCH=R 
GO  TO  4220 

4130 

WAVEFNC I+2)=XPLUS 

4132 

IFIR-RMATCH+.001)  4210,4220,422  0 

4210 

CONTINUE 

4220 

REFL=R*XPLUSP/XPLUS 

WMATCH=XPLUS 

IF(EN)  4211,4216,4216 

4211 

XPLUS=U 

WAVEFN(NP0INT)=1. 

GO  TO  (4212, 4213, 4213), MODE 

4212 

XPLUSP=-SQRTF(VPLUS(NMAX)-EN) 
GO  TO  4214 

4213 

XPLUSP=FLOUT/RMAX 

4214 

H=-H 

R=RMAX 

IN0UT=-1 

DO  5210  J=1,10Q0,KSTEP 

I=NINTVL-J 

CALL  KRR1 

R  =  R  +  H 

WAVEFN( I+1)=XPLUS 

IF(R-RMATCH-.OOl)  4215,4215,5210 

5210 

CONTINUE 

4215 

H=-H 

GO  TO  (5215, 4230), MPRINT 

5215 

1  =  1  +  1 

DO  5216  J=I,NP0INT 

5216 

WAVEFNt J)=WAVEFN( J ) *WMATCH/XPLUS 

4216 

GO  TO  (4221, 4230), MPRINT 

4221 

YINTG=0. 

XINT1=0. 

XINT2=RMAX 

DO  4225  N=1,NP0INT 

4225 

XINTGD(N)=WAVEFN(N)*»2 

CALL  SIMPSN 

DO  4228   N=1,NP0INT 

4228 

WAVEFN(N)=WAVEFN(N)/SQRTF(YINTG) 

4230 

RETURN 

END 
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*  LABEL 
CKRR 

•  LIST 
SUBROUTINE  KRR1 

DIMENSION  WAVE FN ( 20 1),VPLUS( 399), XINTGD (201) 

COMMON  WAVEFN,RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE,NOMAX,NINTVL, 
1L,ELL,K,JDBLD,  JPRINT, MPRINT.NPQINT, EN.FLOUT, HI,R,RR,VPLUS,V,  NMAX, 
2F, RMATCH, I, XPLUS, XPLUSP, FRR, INOUT , K STEP , REFL , XINTGD , YINTG, XI  NT  1 ,H , 
3XINT2,X 
GO  TO  (56,40) ,JPRINT 
40  WRITE  OUTPUT  TAPE  6 ,45 ,H , R,XPLUS, XPLUSP , AI , EN , V, RMATCH 
45  F0RMATL3H  H=,F8.5,3H  R=,F8.5,7H  XPLUS=,  E  12.  5  ,  8H  XPLUSP=  ,E  12.  5  ,4H  A 

1I=,F8.5,4H    EN=,F10.4,3H    V=,F9.4,8H    RMATCH=, F6. 3 ) 
56    RR=R 

INDEX=2*I-IN0UT 

V=VPLUS( INDEX) 

X=XPLUS 

CALL    FR 

AI=FRR*H**2/2.0 

RR=R+H/2.0 

X=XPLUS+    XPLUSP*H/2.0+AI/4.0 

INDEX=2*I+INQUT*(KSTEP-1) 
V=VPLUS(INDEX) 
CALL    FR 

AH    =FRR*(H*»2)/2oO 
RR=R+H 

X=XPLUS+XPLUSP*H+AI I 
INDEX=2*I+IN0UT*(2*KSTEP-1) 
V=VPLUS( INDEX) 
CALL    FR 

AI  II=FRR*H**2/2.0 

XPLUS=XPLUS+H*XPLUSP  +  (AI+2.0*AI  I  )/3.0 
XPLUSP  =XPLUSP  +(AI  +4.0*AII+AIII)/(3.0*H) 
GO  TO  (4125, 190), JPRINT 
190  IF(A8SF(RR-RMATCH)-.001)  195,4125,412  5 

195  WRITE  OUTPUT  TAPE  6 ,45, H, RR ,X PLUS, XPLUSP , AI , EN , V, RMATCH 
4125  CONTINUE 
RETURN 
END 
'*      LABEL 
CFR 

*  LIST 
SUBROUTINE  FR 

DIMENSION  WAV EFN ( 201 ),VPLUS( 399), XINTGD (201) 

COMMON  WAVEFN,RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE,NOMAX,NINTVL, 
1L,ELL,K,JDBLD,-JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR,VPLUS,V,NMAX, 
2 F, RMATCH, I, XPLUS, XPLUSP, FRR, I NOUT,K STEP, REFL, XINTGD, YINTG, XI  NT 1,H, 
3XINT2,X 

FRR=-.04826*F  *(EN     -20.72 1*ELL* ( ELL+1.0) /( F  *RR**2)-V)*X 

RETURN 

END 

*  LABEL 
CSIMPSN 

*  LIST 
SUBROUTINE  SIMPSN 

DIMENSION  W AVEFN ( 20 1),VPLUS( 399), XINTGD (201 ) 

COMMON  WAVEFN,RADIUS,A,VC,.ALPHA,AMASS,RMAX,EPS,MODE,NOMAX,NINTVL, 
lL,ELL,K,JDBLD,.JPRINT,MPRINT,NPOINT,EN,FLaUT,HI,R,RR,VPLUS,V,NMAX, 
2 F, RMATCH, I , XPLUS , XPLUSP ,FRR, INOUT, K STEP , REFL , XINTGD , YINTG, XI  NT  1 ,H, 
3XINT2,X 
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DELINT=(XINT2-XINT1)/FL0ATF(NINTVL) 

NPaiNT=NINTVL+l 

YINTG3=XINTGD(1) 

NPT=l 

DO    20    NPT  =  2,NINTVU2 

YINTG3=yiNTG3+XINTGD{NPT)*4,0 
20    CONTINUE 

DO    30    NPT=3,NINTVL,2 

YINTG3=YINTG3+XINTGD(NPT)*2.0 
30  CONTINUE 

YINTG3=YINTG3+XINTGQ(NP0INT) 

YINTG=YINTG+YINTG3*DELINT/3.0 

RETURN 

END 


20 


Appendix   4.      Fortran   IV   listing. 

A  block  of   COMMON  called  MEMEIG   is   common   to  MAINEV,    EIGENV  and    to 
the    subroutines   called  by  EIGENV.      It  must  appear   in   the   calling   program. 


$IBFTC    MAINEV       LIST,M94,XR7 

C  MAIN    PROGRAM    TO    TEST    EIGENV    SINGLE    PARTICLE    EIGENVALUE    AND    WAVE 

C  FUNCTION    PROGRAM 

DIMENSION    RPRINT(5) , WAVE ( 5 ) , ENRGY ( 1C0),FL( IOC) 

COMMON/MEMEIG/WAVEFN(201),RA0IUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE, 
1N0MAX,NINTVL,L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR, 
2H,VPLUS(399) , V ,-NKAX , F„RMATCH, I ,  XPLUS , XPLUSP t FRR, INOUT, KSTEP, 

3REFL,XINTGD(201) , YINTG, XINT1, XI NT2 , X 
5    WRITE(6,7) 
7    FORMAT! 1H1 ) 
READ(5,10) 
10    FORMAT! 72H 

1  ) 

WRITE(6,10) 
WRITE(6,4501) 
4501    F0RMAK70H  RADIUS  A  VC       ALPHA(VSO)  AMASS 

1RMAX  EPS    ) 

RE AD (5,1)     RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS 
WRITE(6,l)RADIUS,A,VC,ALPHA,AMASS,RMAX,fcPS 
1    FORMAT(7F10.5) 
WRITE(6,4503) 

4503  FORMAT! 70H       ENERGIES  MODE  NOMAX  NINTVL  L  J 
1DBLD            JPRINT     ) 

READ    (5,20)    KZ, MODE, NOMAX  , -NINTVL , L , JD8LD, JPRINT 
WRITE (6, 20)     KZ,MODE,NOMAX,NINTVL,L, J D8LD, JPRINT 
20    FORMAT(7I10) 

NP0INT=NINTVL+1 

GO    TO    (4081, 4082, 4081), MODE 

4081  WRITE(6,4504) 

4504  F0RMAT(15H    TRIAL    ENERGIES     ) 
READ(5,1)  (ENRGY(KY),KY=1,KZ) 
WRITE(6,1)  (ENRGY(KY) ,KY=1,KZ) 
GO    TU    4085 

4082  WRIT£(6,4083) 

4083  FORMAT! 25H   ENERGY        FLOUT,  RMAX  ) 
READ(5,1)    ( ENRGY (KY),FL(KY),KY=1,KZ) 
WRITE (6,1)     (ENRGY(KY),FL(KY) ,KY=1,KZ ) 

4085  DO  100  KY=1,KZ 

EN=ENRGY(KY) 

FLOUT=FL(KY) 

WRITE (6, 6 180) 
6180  FORMAT!  68H  ENERGY     DIFFERENCE      FL  INSIDE       FL  OUT 

1SIDE    RMATCH  ) 

CALL  EIGENV 

IF(ABS     (EN)-l.E-lO)     100,4088,4088 

4088  WRITE(6,4089)  EN 

4089  FORMAT  ( 34H   WAVE  FUNCTION  VS  RADIUS,  ENERGY=, F 10. 5 , 4H  MEV) 
R=-HI 

INMAX=5 

DO  200  N=1,NPCINT,5 

DO  150  IN=1,INMAX 

NN=N+IN-1 
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WAVE(IN)=WAVEFN(NN) 

R=R+HI 

RPRINT( IN)=R 

IF (NN-NPO I  NT) 150,150,15  5 
150  CONTINUE 

GO  TO  159 
155  INMAX=IN-1 

159  WRITc(6,160)  ( RPRINT ( IN ) , WAVE ( I N) , I N=l, INMAX ) 

160  F0RMAT(5(F9.3,E15.6)  ) 
200  CONTINUE 

GO  TO  (100,170), JPRINT 
170  WRITE(6,180)  YINTG 
180  F0RMAT(7H  YINTG= , F 12, 5 ) 
100  CONTINUE 
GO  TO  5 
END 
SIBFTC  EV       LIST,M94,XR7 
SUBROUTINE  EIGENV 
DIMENSION  Y(ll) 

C0MM0N/MEMEIG/WAVEFN(201),.RADIUS,A,VC,ALPHA,AMASS,RMAX,EPS,MODE, 
1N0MAX,NINTVL,L,.ELL,K,JDBLD,  JPRI NT, MPR INT ,NP0INT , EN, FLOUT,HI ,R,RR , 
2H,VPLUS(399)  ,  V  ,.NMAX  ,  F*RMATCH,  I,      XPLUS  ,  XPLUSP,  FRR,  INOUT,KSTEP, 
3REFL,XINTGD(201) , Y I NTG, X  I NT1 , X INT2, X 
K=L+1 

Hl=RMAX/FLOAT(NINTVL) 
NP0INT=NINTVL+1 
IF(ALPHA)  1,2,2 

1  ALPHA=ALPHA*181.1/VC 

2  CALL  VR3 
KSTEP=1 
H=HI 
GO  TO  (15,  15,6130), MODE 

15  IF(EN)  180,180,170 
170  RMATCH=RMAX 

GO  TO  190 
130  RMATCH=(5.»FL0AT(NINTVL/15+1) )*HI 
190  MPRINT=2 

IF(NOMAX-l)  25,30,25 
25  KSTEP=5 

H=5.»HI 
30  CALL  INTEG3 

N0  =  1 
203  ENERGY=EN 

IF(EN)  213,213,220 
213  FLOMT=R*XPLUSP/XPLUS 

DIF=REFL-FL0MT 

GO  TO  6001 
220  PM0M=SQRT(.04826*F»EN) 

DIF=REFL-FLOUT 

FLOMT=FLOUT 

GO  TO  (2201, 6001, 6130), MODE 
2201  XX=R*PMOM 

KPLUS1=K+1 

IFUX-1.)  300,300,400 

300  NPASS=1 

301  SUM=1. 
TERM=1. 
DO  310  N=l,20 

TERM=-TERM*0.5#XX**2/( FLOAT (N)*FL0AT(2«N-1-2*L) ) 
SUM=SUM+TERM 
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IF(ABS(TERM)-l.E-8)  320,310,310 
310  CONTINUE 
320  IF(K-l)  325,335,325 
325  00  330  IK=2,K 

AK=IK 
330  SUM=SUM»(2.*AK-3. ) 
335  GO  TO  (340,350),NPASS 
340  YL=-SUM/(XX**K) 

K=K+1 

L  =  L+1 

NPASS=2 

GO  TO  301 
350  YLPLUS=-SUM/(XX**K) 

K=K-1 

L  =  L-1 

GO  TO  500 
400  Y(1)=-C0S(XX)/XX 

Y(2)=-COS(XX)/XX**2-SIN(XX)/XX 

DO  410  IK=3,KPLUS1 

AK  =  IK 
,  410  Y(  IK)=(2.*AK-3. )*Y( IK-1 ) /XX-Y ( IK-2 ) 
,      YL=Y(K) 

YLPLUS=Y(KPLUS1) 
500  FL0MT=FL0AT(L)+1.-XX*YLPLUS/YL 

DIF=REFL-FL0MT 
5001  IF(NO-NOMAX)  6004,6100,6100 

5004  IF(N0-2)  6005,6080,6090 

5005  EN=0.99*EN 
GO  TO  6100 

5080  EN     =(ENERG1*DIF-ENERGY*DIFM1)/(DIF-DIFM1) 

GO  TO  6100 
5090  EN     =( <iENERG2*DIFMI-ENERGl*DIFM2)*DIF/(DIFMl-DIFM2))-((ENERGl 
1*DIF-ENERGY*DIFM1)*DIFM2/(DIF-DIFM1)) ) / ( DIF-DI FM2 ) 

5100  IF(ABS((EN     /ENERGY)-1. )-,2 )  6110,6110,6101 

5101  EN=ENERGY+ABS( ENERGY ) *S IGN< 0. 1, (EN-ENERGY) ) 
5110  WRITE  (6,6160)  ENERGY, DIF ,REFL, FLOMT , R 
6160  F0RMAT(F15.5,3E15.7,F8.3) 

IF  (NO-NOMAX)  6114,6140,6140 
5114  IF(ABS(OIF     )-EPS  )  6115^6115,6120 

6115  IF(KSTEP-l)  105,6116,105 

6116  EN=ENERGY 
GO  TO  6130 

6120  IF(KSTEP-l)  100,110,100 
100  IF(ABS(DIF     )-.l)  105,110,110 
105  KSTEP=1 

H=HI 
110  CALL  INTEG3 

DIFM2=DIFM1 

DIFM1=0IF 

ENERG2=ENERG1 

ENERG1=ENERGY 

N0=N0+1 

GO  TO  203 
5130  MPRINT=1 

CALL  INTEG3 

GO  TO  6150 
5140  EN=0. 

DO  6145  N=1,NP0INT 
6145  WAVEFN(N)=0. 
6150  CONTINUE 
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RETURN 
END 
$IBFTC  VR       LIST>M94,XR7 
SUBROUTINE  VR3 

COMMON/ MEME I G/WAVEFN  (  20 1),. RADIUS,  A,  VC,  ALPHA,  AMASS,  RMAX,  EPS, MODE, 
1N0MAX,NINTVL,L,ELL,K,JDBLD,JPRINT,MPRINT,NP0INT,EN,FL0UT,HI,R,RR, 
2H,VPLUS(399>  ,  V  ,.NKAX  ,  FrRMATCH,  I  ,      XPLUS, XPLUSP, FRR , IN0UT,KSTEP, 
3REFL,XINTGD(201) , YINTG, XINT1 ,XINT2, X 
R=HI 
H=HI/2. 

GO  TO  (2090,2000)»JPRINT 
2000  WRITEI6,2091) 

2091  FORMAT! 25H  V  R  ) 

2090  C=  .0110270/A 

NMAX=2*NINTVL-1 

DO  2210  I=1,NMAX 

REXP=EXP  UR-RADIUS)/A) 

VCC=VC/(1.0+REXP) 

VS=ALPHA*VC*C*REXP/( ( ( 1 .Q+REXP) »»2 ) *R ) 

AK=K 

IF(JDBLD-2*L)  22,22,21 

21  VPLUSU )=VCC+VS»(AK-1.0) 
GO  TO  25 

22  VPLUS  ( I)=VCC-VS*AK 
25  GO  TO  (40,28),JPRINT 

28  WRITE<6,30)   VPLUSU)*R 
30  FORMATl  E15.7*F7.3) 
40  R=R+H 
2210  CONTINUE 
RETURN 
END 
$IBFTC  INTEG    LIST,M94,XR7 
SUBROUTINE  INTEG3 

COMMON/MEME IG/ WAVEFN( 201), RADIUS, A, VC , ALPHA, AMASS, RMAX, EPS, MODE , 
lNCMAX,NINTVL,L,ELL,K,JDBLD,JPRINTfMPRINT,NPOINT,EN,FLOUT,HI,R,RR, 
2H,VPLUS(399)  ,V„NMAX,F,.RMATCH,  I,      XPLUS, XPLUSP, FRR, IN0UT,KSTEP , 
3REFL,XINTGD(201) , Y INTG, X INT1 , XINT2, X 
GO  TO  (1,15),MPRINT 
1  GO  TO  (5, 5,15), MODE 
5  IF ( EN)  5215,4221,4221 
15  WAVEFN(1)=0. 
AK  =  K 

F=AMASS/( AMASS+ 1.00898) 
R=H 
RR  =  R 

V=VPLUS(2»KSTEP-1) 
XPLUS=R«*K 
X=XPLUS 
ELL=0. 
CALL  FR 
ELL=K-1 

XPLUS=XPLUS+FRR*R«»2/((ELL+2. )*(ELL  +  3.)  ) 
WAVEFN(2)=XPLUS 

XPLUSP    =AK*(R  **(K-D)  +FRR*R/ (  ELL  +  2.  ) 
GO  TO  (4090, 4090,4080), MODE 
4080  N0MAX=2 
4090  IN0UT=1 

DO  4210  I=KSTEP,1000,KSTEP 

1=1 

CALL  KRR1 
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R  =  R+H 

4127  IF(NOMAX-l)  4130,4128,4130 

4128  RMATCH=RMAX 

IF  (R-0.2*RADIUS)  4130,4129,4129 

4129  IF(EN-V)      10,10,4132 
10  RMATCH=R 

GO  TO  4220 

4130  WAVEFN(I+2)=XPLUS 

4132  IF(R-RMATCH+.001)  4210,4220,4220 

4210  CONTINUE 

4220  REFL=R*XPLUSP/XPLUS 
WMATCH=XPLUS 

IF(EN)  4211,4216,4216 

4211  XPLUS=1. 
WAVEFN(NP0INT)=1. 

GO  TO  (4212, 4213, 4213), MODE 

4212  XPLUSP  =-SQRT(VPLUS(NMAX)-EN) 
GO  TO  4214 

4213  XPLUSP=FLOUT/RMAX 

4214  H=-H 
R=RMAX 
IN0UT=-1 

DO  5210  J=1,1000,KSTEP 
I=NINTVL-J 
CALL  KRR1 
R=R+H 

WAVEFN( I+1)=XPLUS 

IF(R-RMATCH-.OOl)  4215,4215,5210 
5210  CONTINUE 

4215  H=-H 

GO  TO  (5215, 4230), MPRINT 

5215  1=1+1 

DO  5216  J=I,NPOINT 

5216  WAVEFN(J)=WAVEFN(J)*WMATCH/XPLUS 

4216  GO  TO  (4221, 4230), MPRINT 

4221  YINTG=0. 
XINT1=0. 
XINT2=RMAX 

DO  4225  N=1,NP0INT 
4225  XINTGD(N)=WAVEFN(N)**2 

CALL  SIMPSN 

DO  4228   N=1,NP0INT 
4228  WAVEFN(N)=WAVEFN(N)/SQRT(YINTG) 

4230  RETURN 

END 

iSIBFTC  KRR      LIST 

SUBROUTINE  KRR1  _  „u4M  cnc    unr.c 

COMMON/MEMEIG/WAVEFNf 201), RADIUS, A, VC, ALPHA, AMASS, RMAX, EPS, MODE, 

1N0MAX,NINTVL,L^ELL,K,JDBLD,JPRI NT, MPRINT, NPOINT, EN, FLOUT,HI,R,RR, 

2H,VPLUS(399),V,NMAX,FrRMATCH,I,      XPLUS , XPLUSP, FRR , INOUT,KSTEP, 

3REFL,XINTGD(201) , Y INTG, XINT1, X  I NT2, X 
GO  TO  (56,40) ,JPRINT 
40  WRITE(6,45)    H, R, XPLUS , XPLUSP, AI , EN, V, RMATCH 
45  F0RMAT(3H  H=,F8.5,3H  R=,F9.5,7H  XPLLS=, E12. 5, 8H  XPLUSP= ,E12. 5,4H  A 

1I=,F9.5,4H  EN=,F10.4,3H  V=,F9.4,8H  RMATCH=, F6.3 ) 

56  RR=R 

INDEX=2*I-IN0UT 
V=VPLUS( INDEX) 
X=XPLUS 
CALL  FR 
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AI=FRR*H**2/2.0 

RR=R+H/2,0 

X=XPLUS+    XPLUSP»H/2.0+AI/4.0 

INDEX=2*I+IN0UT*(KSTEP-1) 

V=VPLUS(INDEX) 

CALL    FR 

AH    =FRR»(H**2)/2.0 
RR=R+H 

X  =  XPLUS+XPLUSP*H+AI  I 

1NDEX=2*I+IN0UT*(2*KSTEP-1) 

V=VPLUS(INDEX) 

CALL    FR 

AI II=FRR*H**2/2.0 
XPLUS=XPLUS+H»XPLUSP+(AI+2.0*AII)/3a0 

G0LTOP<:i2SUiqPn;(AI    **-0-*II*Ai«I?lS.O«l 

ku    TO    (4125,190) ,JPRINT 

190    IF(A8S(RR-RMATCH)-.00D     195,4125,4125 

RETURN 
END 
$IBFTC    FR  LIST*M94,XR7 

SUBROUTINE    FR 
C0MM0N/MEMEIG/WAVEFN(2Ol >    RAnnic    a    ur    *.„,. 

3REFL,XINTGDI201),YINTG?«;n:«NT2.x  '  XPLUSP'FRR' '"OUT,  KSTEP  . 

RE?WN°4826*F    *<EN  -^WJl.ELLMELL^.o./IF    .RR»2,-v,.X 

END 

SIBFTC  SIMPSN   LIST,M94,XR7 
SUBROUTINE  SIMPSN 

DEL INT=(XINT2-XINT1) /FLOAT  (NINTVL) 

YlNrG3=XINTGD(l) 

NPT=1 

DO  20  NPT=2, NINTVL, 2 

YINTG3=YINTG3+XINTG0(NPT)#4.0 
20  CONTINUE 

DO  30  NPT=3, NINTVL, 2 

YINTG3=YINTG3+XINTGD(NPT)*2.0 
30  CONTINUE 

YINTG3=YINTG3+XINTGD(NP0INT) 

YINTG=YINTG+YINTG3*DELINT/3.C 
RETURN 

END 


-EPS, MODE, 
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Figure  Captions 

Results  of  the  "signature  of  the  well"  run  for  the  £=0   well 
shown  in  Figure  3.   Results  are  given  for  f  (inside),  f  (out- 

Xj  Hi 

side),  and  for  the  difference  (DIF)  between  them.   Note  that 

the  energy  at  which  the  difference  is  zero  is  close  to  the 

energy  of  the  converged  eigenvalue,  indicated  by  the  short 

arrow. 

Improvement  of  the  wave  function  match  during  the  automatic 

search  procedure.   Note  that  for  proper  normalization  of  the 

bound  state  wave  function  that  RMAX  must  be  sufficiently  large 

2 
that  there  is  negligible  contribution  to  m   dr  beyond  RMAX. 

Single  particle  states  and  the  lsj,  2si,  and  ld-5  wave  functions 

8      2  2 

in  a  Woods-Saxon  potential  well.   The  parameters  of  the  well 

are:  V  =  -50  MeV,  RADIUS  =3.15  fermis,  a  =  0.650  fermis,  and 
c 

the  spin-orbit  term  is  20  times  the  Thomas  term  for  a  nucleon 
(ALPHA  =  20).   The  shape  of  the  well  shown  is  for  4=0.   During 
the  "signature  of  the  well"  preliminary  run,  the  matching 
radius  is  given  by  the  heavy  black  line  -  the  edge  of  the  well 
for  bound  states,  and  the  maximum  radius  (8  fermis  in  this 
example)  for  continuum  states. 
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